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Abstract 



We consider the aggregation equation pt—V-{pVK * p) = in M"', where the interaction 
potential K incorporates short-range Newtonian repulsion and long-range power-law 
attraction. We study the global well-posedness of solutions and investigate analytically 
and numerically the equilibrium solutions. We show that there exist unique equilibria 
supported on a ball of M". By using the method of moving planes we prove that 
such equilibria are radially symmetric and monotone in the radial coordinate. We 
perform asymptotic studies for the limiting cases when the exponent of the power-law 
attraction approaches infinity and a Newtonian singularity, respectively. Numerical 
simulations suggest that equilibria studied here are global attractors for the dynamics 
of the aggregation model. 

Keywords: swarm equilibria, biological aggregations, Newtonian potential, global at- 
tractors 



has attracted a great amount of interest in recent years. The equation appears in various 
contexts related to mathematical models for biological aggregations, where p represents the 
density of the aggregation and K is the social interaction potential. The asterisk * denotes 
convolution. We refer to [321 [M] for an extensive background and literature review on math- 
ematical models of social aggregations and in particular, for a thorough discussion on the 
relevance of equation ([T]) for modelling swarming behaviours. The equation also arises in 
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1 Introduction 



The multidimensional integro-differential equation. 



Pt-V- {pVK * p) = 0, 
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a number of other applications such as granular media [351 IlZj, self-assembly of nanopar- 
ticles p4l [25] . Ginzburg-Landau vortices [T6l [T4l 131] and molecular dynamics simulations 
of matter [23] . In this work however we are primarily interested in biological applications, 
where equation ([T]) is used to model social aggregations such as insect swarms, fish schools, 
bacterial colonies, etc [32] . 

Regarded as a model for biological aggregations, equation ([1]) incorporates inter-individual 
social interactions such as long-range attraction and short-range repulsion, through the ag- 
gregation potential K. The properties of the potential (symmetry, regularity, monotonicity, 
etc) are essential in studying issues such as the well-posedness [HI El [8] or the long-time 
behaviour [TOl [28] of solutions to model equation ([1]). In particular, a large component 
of the research on this model dealt with attractive potentials K which lead to solutions 
that blow-up (in finite or infinite time) by mass concentration, into one or several Dirac 
distributions [Mlil[25]. 

It is essential however for an aggregation model to be able to capture solutions with 
biologically relevant features. As pointed out by Mogilner and Keshet in their seminal 
work [32] on the class of models discussed here, such desired characteristics include: finite 
densities, sharp boundaries, relatively constant internal population and long lifetimes. The 
difficulty in finding such solutions to model ([T]) has been indicated as a "challenge" in previous 
literature [331 [29], ^i-nd in fact there is only a handful of works that address this issue. Topaz 
and collaborators [281 Ej derived explicit swarm equilibria that arise in the one- dimensional 
model with Morse-type potentials (in the form of decaying exponentials), but their explicit 
calculations do not extend to higher dimensions. Other works illustrate asymptotic vortex 
states in 2D [33] and clumps (aggregations with compact support) in a nonlocal model that 
includes density- dependent diffusion [M]. 

A recent publication of the authors [2U] considered an interaction potential K for which 
equilibria of the aggregation model ([1]) have the desired characteristics indicated above. 
More specifically, the kernel investigated in [20] has a repulsion component in the form of 
the Newtonian potential and attraction given by a power 



Here, 0(x) is the free-space Green's function of the negative Laplace operator —A: 



and g is a real exponent, q >2. In ([3]), n is the number of space dimensions and Un denotes 
the volume of the unit ball in R"". 

We summarize briefly some of the results from [20] that are relevant to the present article. 
For q = 2, the equilibrium density of ([I])-([2l) is uniform inside a ball of M" and zero outside 
it. In this case, the method of characteristics was used to solve explicitly the dynamics 

^See Section [5] for a discussion on how the potential can be modified to avoid the biologically unrealistic 
growth of attraction with distance when q > 0. 



K{x) = 0(x) + 



(2) 
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corresponding to radially symmetric initial conditions in any dimension. This showed the 
global stability within the class of radially symmetric solutions of the constant steady state. 
The explicit calculations did not extend to general exponent q > 2, but the existence of 
a unique radially symmetric equilibrium of compact support was shown, after casting the 
equilibrium problem as an eigenvalue problem for an integral operator and applying the 
Krein-Rutman theorem. Some explicit calculations of the equilibria could be performed 
however for the special subcase when q is even. In addition to these studies on equilibria, 
the global well-posedness of solutions to ([I])-(l2]) (with q > 2) was shown by borrowing 
techniques used in the analysis of incompressible Euler equations [30] . 

The main purpose of the work from [20] was to design attractive-repulsive potentials 
that yield equilibrium states of finite densities and compact support. In this respect, the 
attraction component ^\x\'^ of the potential was specifically designed to counter-balance 
the singular Newtonian repulsion. Newtonian (attractive) potentials for model ([1]) were 
considered in [T^ IHTj in the context of vortex motions in two-dimensional superfluids. The 
main concern of these works was the well-posedness of solutions, in particular concentration 
and singularity formation in measure-valued solutions. Very recently, Newtonian potentials 
were also considered in aggregation models [3 [7]. In [7], the authors study patch solutions 
and they consider separately the case of an attractive Newtonian potential (with finite time 
concentration) and of a repulsive Newtonian potential (with spreading to a circular/ spherical 
aggregation patch). 

The purpose of the present research is (i) to extend the interaction potential (I2])-(|3]) 
studied in [20] to allow for more general attractive forced {q > 2 — n) and (ii) to investigate 
analytically and numerically the properties of the equilibria to the aggregation model ([1])- 
for g > 2 — n. Remarkably, the intricate balance between the power-law attraction 
and the singular repulsion provides the model with a very interesting and at the same time 
biologically relevant set of steady states. For all values of g G (2 — n, oo), the aggregation 
model has a unique steady state supported in a ball. This steady state is radial and monotone 
in the radial coordinate. More specifically, the equilibria are decreasing about the origin for 
2 — n < q < 2 and increasing for q > 2, while q = 2 corresponds to a constant equilibrium 
density. Figure [T] shows the equilibrium solutions in three dimensions for various values of q; 
all shown equilibria have mass 1. The limits g — )■ oo and g \ 2 — n, that is, when attraction 
becomes infinitely strong (at large distances) or as singular as the (Newtonian) repulsion, are 
particularly interesting. As g — ?■ oo, the radii of the equilibria approach a constant, but the 
qualitative features change dramatically, as mass aggregates toward the edge of the swarm, 
leaving an increasingly void region in the centre — this effect can be observed in Figure [T] 
(g = 20,40,80). As g \ 2 — n, the radii of equilibria approach and mass concentrates 
at the origin - see Figure [1] (g = 1.5,1,0.5). Numerical simulations suggest that all these 
equilibria are global attractors for the dynamics of ([I])-([2]), which motivates and gives strong 
grounds to the studies of the present work. 

There are a few very recent studies of equilibria of ([T]) with attractive-repulsive potentials 
in power-law form that closely relate to our work. In [T], the authors study the stability of 
spherical shell equilibria for potentials in the form K{x) = — where 2 — n<p<q 
(short range repulsion and long range attraction). While the attraction component, \x\'^/q, 

^In the special case g = we take K{x) = (/)(x) + In |a;|. 
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Figure 1: Radially symmetric equilibria of (JT])-® in three dimensions, for various values 
of g.The equilibria are monotone in the radial coordinate: decreasing about the origin for 
2 — n < g < 2, increasing for q > 2, and constant for q = 2. As g — >■ oo, the radii of the 
equilibria approach a constant, and mass aggregates toward the edge of the swarm, leaving 
an increasingly void region in the centre. As g \ 2 — n, the radii of equilibria approach 
and mass concentrates at the origin. Numerics suggests that all these equilibria are global 
attractors for the dynamics of ([I])- ([2]). 

is the same as in ([2]), the singularity of the repulsion term is "better" than Newtonian 
{p > 2 — n). Shell steady states are shown to exist and be locally stable under certain 
conditions on the exponents p and q. The methods from [1] do not apply to Newtonian 
singularities. Other recent works that involve model ([1]) with power-law potentials focus on 
pattern formation and linear stability analysis of spherical shells [271 EZ] • 

The results of the paper are as follows. Well-posedness of solutions to the aggregation 
model ([I])-© (with q > 2 — n) is studied in Section[2]by using analogies with the incompress- 
ible fluid flow equations [301 EDI E]- For all values of g G {2 — n, oo), we show in Section [XT] 
that there exist unique equilibria supported on a ball of M". In Section 13.21 we employ the 
method of moving planes [21] to prove that such equilibria are radially symmetric and mono- 
tone in the radial coordinate. In Section |4] we performed careful asymptotic and numerical 
investigations of the equilibria. Our studies address two issues. The first one is the be- 
haviour of equilibria as g — ?■ oo and g — ?■ 2 — n. As expected, the two limiting cases give very 
different asymptotic behaviours. The second issue addressed in Section H] is the stability of 
the equilibrium solutions. The results regarding stability are preliminary and entirely based 
on numerical observations. Based on all numerical experiments we performed, we conjecture 
that the equilibria studied in this paper are global attractors for solutions to (II])-®. 
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2 A priori bounds and well-posedness of solutions 

We start by pointing out that the aggregation model ([I])-(l2l) has two important conservation 
properties. Denote the initial density by po- 

p{x,0) = po{x), xeW. 

The aggregation model ([I])-® satisfies: 

(i) Conservation of mass: 

J p{x, t)dx = M, for all t > 0, (4) 

where the constant M denotes the initial mass M = J pQ{x)dx. 

(ii) Conservation of centre of mass: 

J xp{x,t)dx = 0, for all t > 0, (5) 

where we assume, without loss of generality, that the centre of mass of the initial density is 
at the origin: f xpo{x)dx = 0. 

Both properties follow directly from ([1]). Property (ii) uses the radial symmetry of the 
potential. The two conservation properties will be used frequently in this article. 

By introducing the notation: 

fix) = -VK{x), 



we write the aggregation model as 

Pi + V ■ (pv) = 0, (6a) 
v = f* p, (6b) 

\nujn \x\"- ^ J \x\ 

This work makes extensive use of the Lagrangian formulation of the aggregation model 
where dynamics is tracked along the characteristic curves, defined by: 

^X{a,t) = v{X{a,t),t), X{a,0) = a, (7) 
with velocity v defined by (]6b|) and ( pc|) . 

2.1 A priori bounds on density 

Expand V ■ (pv) = v ■ Vf + pV ■ v and write the evolution of the density p{X{a, t),t) along 
characteristics: 

^{X{a,t),t) = -pV ■ v{X{a,t),t). (8) 



Calculation of V ■ w from ( pojl and (l6c|) yields: 

V-v = p-{n + q-2)l \x - y\'^^'^ p{y)dy. (9) 

The continuity equation ( l6all expresses the fact that 

p(X(a,t),t)J(a,t) = po(a), (10) 

where 

J{a,t) = det V„X 

is the Jacobian of the particle map a — )■ X{a, t). 
Define the maximum of the density 

Pmax(t) = SUpp(x,t). 

X 

We show that provided pmax is bounded initially, it remains bounded above, uniformly in 
time. This was shown in |20] for q > 2 and here we extend the results to include 2— n < g < 2. 

Let q E (2 — n, 2). From the characteristic equation ([8]) for p and the expression ([9]) for 
V • w, we have, along particle trajectories: 

^ = -p' + in + q-2)p [ \x-y\'^''piy)dy (11) 
at 



= -p'^ + {n + q- 2)p 
< -p^ + {n + q-2)p 



|>r* 



x-y\'' piy)dy 



/ Ix-yl" '^p{y)dy+ / 

.J\x—y\<rt J\x—y 

Pmax / \x- y\'^~^dy + rl'^ / p{y)dy 

J\x—y\<rt J\x—y\>rt 



(12) 



where we used q < 2 and r^, > will be chosen conveniently later. Use the following estimates 
on the integrals in the right-hand-side of flT^ : 



/ 



x-yr'dy= r^^-^ / p{y)dy<M, 



and choose to be: 

r* = (Prnax)-'/". (13) 

We find 



< -p2 + Cpp(^^^ , 



dp 
with 

C = + M{n + g - 2), 
resulting in the following inequality 



^ < C^pL^aic^"^^/" - P^a.- (14) 
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As g > 2 — n, the damping dominates the growth term in the right-hand-side of (fT4|) : 
{n + 2 — q)/n < 2. The right-hand-side of ( fT4l) becomes negative when pmax > C^"'^''~'^^^"' 
and hence, regardless of the size of the support, the maximum density is bounded uniformly 
in time, provided it is initially bounded. 

Remark. For g > 2, it was shown in [20] that the density has compact support uniformly 
in time, provided the initial density po has compact support. This property was used to 
conclude global wellposedness of solutions. In the present study, where 2 — ?7,<g<2, we 
could not show the compact support of solutions when 2 — n < q < 1, but we managed to 
circumvent this by using the uniform L^-bound of p. 

We present briefly the argument that shows uniform compact support for 1 < g < 2. 
The density is transported along characteristics (see equation ffTOl) ). so it is enough to show 
that the trajectories X{a,t) that carry non-zero densities remain within some compact set. 
Calculate using ( l6b|) and ( l6cl) : 

f X ' (x — f 

x-v{x,t)= ■ -p{y,t)dy- x ■ {x - y)\x - y\''^^p{y,t)dy. (15) 

Jm" nunix — y\ 

Define the maximum radius of support R{t) as 

R(t) = max |X(a,t)|, 

a:po{a)^0 

and evaluate ( fTSll at x = X{a,t) on the boundary of the support, i.e., |X(a,t)| = R(t) > 
\X{f3,t)\ for any (3 such that po{f3) > 0. The left-hand- side of (fT5|) becomes 

X(a,t) ■ ^X(a,t) = R^. 
^ ' ' dt ^ ' ' dt 



We estimate the first term in the right-hand-side of (1151) as follows: 

x-(x — ?/),,, R f 1 /N, R f 1 , \ 1 
■ r^p{y,t)dy< / n^p{y)dy+ / —p[y)dy 

M 

< PniaxR H R- 

nUJn 

For the second term, use |x — y| < 2R and x ■ (x — > 0, for |x| = i? on the boundary of 
the support and y in the support of p, to find 

X ■ (x - y) |x - i)dy < -{2Ry~'^ / X ■ (x - y)p{y)dy 

= -2''-'^MR\ 

where we used conservation of mass (jlj) and centre of mass ^ to go from the second to the 
last line. Using the estimates in (1151) . derive 

dR M „ r, 1 

-1- <Pmax + 2''~^MR'^~\ 

dt noOr, 



^max, where 



Hence, the trajectories that carry non-zero densities will remain inside the disk of radius 

_i 

_ ■ Pmax + M/ (nUJn ) \ " 
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2.2 Existence and uniqueness of solutions 

To study well-posedness of solutions we use a Lagrangian approach and rewrite the aggrega- 
tion equation ([T]) in terms of particle trajectories. Then we regard the model as an ODE on 
a certain Banach space and infer local existence and uniqueness from Picard theorem. The 
setup of the ODE framework is the same as that used in [20] to study the case q > 2 and is 
inspired from the study of well-posedness of solutions to the incompressible Euler equation in 
Lagrangian formulation [30]. Extension to global existence is achieved through an argument 
similar to the well-known Beale-Kato-Majda blow-up criterion for incompressible flows [2]. 

Make the change of variable y = X{fi, t) in the expression (I6b|l for f , with / given by (l6c|l . 
and use (fTOl) to write the characteristic equation ([7]) as 

^Xia.t) = J^iXia.t)) (16a) 
at 

X(a,0) = a, (16b) 

where the map J^{X) is defined by 



nUn \X{a,t) - X{l3,t)\ 

(17) 

System (fT6l) -( |T71) is a reformulation the PDE model ([6]) in terms of particle-trajectory 
equations. Case q > 2 was studied in detail in [20], by analogy with the ODE setup of the 
incompressible Euler equations [30] . 

Following [301 l2D] , we consider the Banach space 



= {X : M" ^ such that < oo}, 

where || ■ is the norm defined by 

= |X(0)| + ||V«X||lco + \VaX\y. {li 

Here, | ■ | is the Holder seminorm 

,^ ^, \VaX(a)-\/r.X(a')\ 
I VqA = sup 

Consider an open subset Ol, of B defined by 



\a — a 



'|7 



Or 



|x e I inf det V„X(a) > l/L and < l| . 



The key ingredients to show local and global well-posedness of solutions are the properties 
of the convolution kernel 

k{x) = (19) 

present in the repulsion component of ( lUcj) . In particular, k is singular, homogeneous of 
degree 1 — n and its gradient P = Vk is homogeneous of degree —n and defines a singular 
integral operator (SIO). The close analogy with incompressible fluid equations comes from 
the fact that a similar kernel appears in the Biot-Savart law [30] . 

The local existence and uniqueness is stated by the following theorem. 
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Theorem 2.1. (local existence and uniqueness) Consider a compactly supported initial 
density po G L°°(M.'^), with |pol7 < oo, for some 7 G (0, 1). Then for any L > 0, there exists 
T{L) >0 and a unique solution X E C^{{-T{L),T{L)); Ol) to with q > 2 - n. 

Proof. Following [201 EHj , we show that the map J-" is bounded and locally Lipschitz contin- 
uous on Ol- Local existence and uniqueness then follows from Picard theorem. Details are 
presented in the Appendix. □ 

We now use a continuation result of solutions to autonomous ODE's on Banach spaces 
(Theorem 4.4 in [30]) to upgrade the result to global existence. Inspecting the set Ol we 
infer that we cease to have a solution at a finite time provided either inf^j det VQX(a) 
becomes or becomes unbounded as t — T*. 

The first scenario is ruled out by the following proposition. 

Proposition 2.2. At any fixed time t < 00, solutions of (IT6l)-( fT71) satisfy 

inf det V„X(a,t) > e"^*, 

a 

where C > is a constant. 

Proof. Use the differential equation that J satisfies, 

d 



to derive 



-J{a,t) = J{a,t)V ■ v{X{a,t),t), 



J(a,t) = exp^y V ■ v{X{a, s), s)ds^ . 



Case q > 2 was considered in [20]. For 2 — n<q<2, one can use and similar calculations 
to those leading to ffT^ - (fn|) . to derive 

2-q 

|V ■ w| < IIpIIl- + C'IIpII^So . 
As ||p||l°° is uniformly bounded in time (see Section I2TT]) . we find 

J{a,t) > exp(-Ct), 

with C > 0. □ 

The second scenario for the break-up of the solution (finite-time blow-up of will 
be treated as in Chapter 4.2 of [30]. This procedure was used in [20] to study global well- 
posedness of solutions to ( fT6l) . ( fTTI) when q > 2. In summary, can be shown to remain 

bounded for all finite times, provided \\p{-, s)ds\\ds < 00, for all t. This is an analogue 
of the Beale-Kato-Majda condition for global existence of solutions to incompressible Euler 
equations [2], [30]. Our argument is adapted from the analysis of the incompressible fluid 
equations presented in Chapter 4 of [30] . 

A first a priori bound is provided by the following proposition. 
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Proposition 2.3. Provided /q || Vw(-, s)||ioo has an a priori bound, || Vq,X(-, t)||ioo and |VqX(-, 
are a priori bounded. 

Proof. The proof for q = 2 was presented in detail in [20] (see Proposition 2.3). Extending 
it to general q > 2 — n does not pose any difficulties and we omit the details. See also 
Proposition 4.3 from [30] for the corresponding result in the context of incompressible Euler 
equations. □ 

The Beale-Kato-Majda condition for global existence of incompressible Euler equations 
is an a priori control on the time-integral of the supremum norm of vorticity. In the context 
of our aggregation model fll6p - (IT7|) . this condition will be replaced by an a priori bound on 
/o IIp(->s)IIl- ds. 

Proposition 2.4. A sufficient condition for ||V?7(-, s)||ioo(is to be a priori bounded is an 
a apriori bound on \\p{-, s)\\Lcx,ds. More specifically, 

«+ 

where C{po) is a constant that depends on the initial density only. 

Proof. Case q > 2 was discussed in [20]. The proof for 2 — n < q < 2 requires some slight 
adaptations from the corresponding result for fluids (see Theorem 4.3 and its proof in |30j). 
The repulsion component of Vt' can be estimated using (Ell). The attraction part has a 
milder singularity and does not break the estimate, hence we have 

II Vi;(-, t) lUoo < c (|p(-, t) l^e^ + ||p(-, t) lUco log(l/e)) + M. 

Set e = |p(-,t)|7^/^ to get 

||V^;(-,t)|U- < ||p(-,t)||L- (log|p(-,t)|^ + c). 



Lemma 4.8 in [30] can be trivially adapted to our context, resulting in the following inequal- 
ity: 

|p(-,i)l7 < IPo^exp ((^ + ^)^ ||Vt;(-,s)||Loods^ . 
By combining the last two inequalities we find 

||V^;(-,t)|Uoc <C(po)||p(-,t)|U- (^1 + ^ ||Vt;(-,s)|Uocrfs^ . 

The desired inequality follows after division by 1 + Jq* || Vw(-, s)||Loo(is and integration with 
respect to t. □ 

Finally, we have all the ingredients to prove global existence of solutions. The result is 
given by the following theorem. 
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Theorem 2.5. (global existence and uniqueness) Consider the trajectory equations (fT6|) . 
(fTTl) with the Banach space setup and notations as above, and a compactly supported initial 
density po G L°°{MP), with \po\^ < oo, for some 7 G (0, 1). Then, for every T, there exists 
L > and a unique solution X G C^{[0,T)] Ol) to ( |T6l) . ( |T7j) with q > 2 — n (a unique 
solution exists globally in time). 

Proof. Case q > 2 was studied in [20] and we focus here on the range 2 — n < q < 2. 
The solution X is in the set Ol provided 

midetVaX{a) > 1/L and < L. (20) 

Using Proposition 12.21 the first condition is satisfied provided we choose L > e*"^. We now 
investigate the second condition in ( l20i) . Start by inspecting the first term in ( fT8l) . |X(0,t)|. 
Integrate to get 

X(0,t) = / v{X{0,s),s)ds. (21) 



The repulsion component of v can be bounded in terms of ||p||lo° using ([SSD- The attraction 
component also has a uniform bound. To show this we distinguish two cases: (i) 2—n < q < 1 
and (ii) 1 < g < 2, and inspect the attraction part of v, i.e., — * p. In case (i), 



Ixl" ^x*p\< [ \—p{y)dy+ [ \—p{y)dy 

J\x^y\<i\x-yr~^ J\x-y\>i\x - y^^'" 

< c\\p\\lo^ + M. 



In case (ii) solutions p are compactly supported (uniformly in time) in a ball of radius R, 
provided the initial density po is — see Remark in Section 12.11 Hence, 



p\ < {2Rf-^M. 



Using fl2T|) we derive 

\X{0,t)\<ci [ \\p{-,s)\\Lo.ds + C2t, (22) 
Jo 

The control of now follows from fl22|) . Propositions 12.31 and 12.41 More precisely, by 

redefining the constants appropriately, one can derive 

||X(-,t)||i,,<Ce^^V^<^^"'°'"''-'^"^°°". (23) 

Given that ||p(-, s)||2,oo is uniformly bounded, we can choose the constant L large enough 
such that ||X(-,t)||i,^ < L, for aU t G [0,T). 

□ 



3 Radially symmetric steady states 

In this section we show that the aggregation model admits a unique radially symmetric 
steady state p supported on a ball of M". For 2 — n < g < 2, we prove that these steady 
states are monotonically decreasing about the origin, while for q > 2 they are increasing 
about it. Case q = 2 corresponds to constant solutions in a ball [20]. We further study these 
equilibria using numerical and asymptotic methods in Section HI 
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3.1 Existence and uniqueness of equilibria supported on a ball 

Suppose that p is a steady state with support the ball -8(0, R) centred at the origin, of radius 
R. The velocity v is zero in -8(0, R), so its divergence also vanishes. Hence, from ([9]) we find 
that p satisfies the following integral equation, 

p{x) -{n + q-2) I \x-y\'i~'^p{y)dy = for x G 5(0, -R), (24) 

and vanishes outside B{0,R). 

Consider the operator Tr given by 

Tnp{x) = {n + q-2) [ \x - y\'^-''p{y)dy. (25) 

Jb{o,r) 

The subscript is used to emphasize the dependence of the integral operator on the radius R. 
As g — 2 > — n, the kernel \x — y\'^~'^ is integrable and Tr defines a linear bounded operator 
from C(-B(0, -R), M) to itself. Equation (!24l) can be cast as an eigenvalue problem, 

TRp=-p, pGC(fi(0,i?),M), 

where solutions p are eigenfunctions corresponding to eigenvalue 1. 

The existence and uniqueness of a steady state supported in -6(0, R) is provided by the 
following theorem. 

Theorem 3.1. For every q > 2 — n and M > 0, there exists a unique radius R (that 
depends on q and n only) and a unique steady state p of the aggregation model that 
is supported on B{0,R), has mass M and is continuous on its support. 

Proof. We use a scaling argument and consider the case -R = 1 first. For q > 2 — n, 
the kernel \x — is integrable and the linear operator Ti : (C(-B(0, 1), R), || ■ ||l°°) 

{C{B{0, 1), M), II ■ ||l°°) is bounded. The operator is also compact. This is a textbook exercise 
in analysis [36], but we include it here for completeness. Case q > 2 presents no difficulty, 
as the kernel \x — ?/|^~^ is continuous. For 2 — n < g < 2, we write Ti as a limit of compact 
operators as follows. Consider a smooth cut-off function h : [0, oo) — t- M, such as h{r) = 
for < r < 1/2, h{r) = 1 for r > 1 and < h{r) < 1 for all r > 0. Define 

Kp{x,y) = h{p\x -y\)\x -yl"-'^, 

and 

TipP{^) = {n + q-2) / Kp{x, y)p{y)dy. 

Jb(o,i) 

The kernels Kp are continuous, hence the operators Tip are compact. Now estimate 

\Tip{x) - Ti(x)| < {n + q-2) [ \x - y\''~^\l - h{p\x - y\)\p{y)dy 

Jb(oa) 



mi) 

< {n + q - 2)\\p\\l^ / \x-y\''~'^dy 

J\x-y\<l/p 



1 

< nujn\\p\ 



pn+q—2 ' 
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As n + g — 2 > 0, ^„^\„2 0, when p — t- oo. The convergence is uniform in x, hence, 

\\Tip - Till ^ 0, p~> oo, 

and Ti is compact as a hmit (in the operator norm) of compact operators. 

We apply the Krein-Rutman theorem [15] to operator Ti in the following setup. Take the 
cone in C{B{0, 1), M) consisting of all non-negative functions. Ti is a linear, strongly positive, 
compact operator that maps the space of continuous functions C{B{0, 1),M) into itself. By 
Krein-Rutman theorem (see Theorem 1.2 in [I5]), there exists a positive eigenfunction pi 
such that Tipi = Api, where A (which depends only on q and n) is the spectral radius of 
Ti. Moreover, the eigenvalue A is simple and there is no other eigenvalue with a positive 
eigenvector. By making the change of variable 

p{x) = pi{x/R) (26) 

in (125|) . we get 

Tnp{x) = R'+'^-^Xpix). 
Now ask that p is an eigenfunction of Tr corresponding to eigenvalue one and find 

R = X-^^^ (27) 

which gives the radius of the support as a function of q and n only. Once a mass M for p 
is set, uniqueness can be inferred from the uniqueness properties of the spectral radius of Ti 
and its associated eigenfunction pi. □ 



3.2 Qualitative properties of equilibria 

Monotonicity and radial symmetry of equilibria. We prove that the equilibria given 
by positive solutions of ( 124|) (and whose existence and uniqueness was established in The- 
orem [STT]) are radially symmetric and monotone in the radial coordinate. More specifically, 
equilibria are monotonically decreasing when 2 — n < q < 2 and increasing for q > 2. To 
prove this result we employ the method of moving planes, a technique introduced by the So- 
viet mathematician Alexandroff in the early 1950's, which became well-known after Gidas, 
Ni and Nirenberg [2T] applied it to study qualitative properties of positive solutions of elliptic 
equations. We refer to [13] for a detailed description of the method and its applications. 

Our use of the moving plane technique is inspired by a novel application of the method 
in the context of integral equations [13]. Consider a steady state p supported on the ball 
B{0,R); p satisfies the integral equation fl2^ in B{0,R) and vanishes outside B{0,R). For 
convenience of calculations, denote 

a = n + q — 2. 

As g > 2 — n, we have a > 0. Two cases will be distinguished from the subsequent analysis: 
(i) 2 — 71 < g < 2 (equivalently, < a < n) with radial equilibria which decrease about 
the origin, and (ii) g > 2 (or a > n) with increasing solutions. Using the new notation, p 
satisfies 

p(^) = Ib{o,r) dy xe B{0, R) ^^^^ 

[0 x^B{0,R). 
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Take /i G M such that —R < /i < 0, and consider the reflection across the plane xi = jj,, 
X H-> = (2/i — xi,X2, ■ ■ ■ ,Xn)- In particular, we have the image O'^ of the origin under 
this map. Define 

Pf,{x) = pix"). 

Using ( 128|) . is given by 

^^(3.) = 1" lBio,R) W^\^Piy) dy xe 5(0^ R) ^^g^ 



x^B{0^',R). 



Define 



= {x e B{0,R) I xi > p}. 

We apply the method of moving planes and compare p{x) and Pfj,{x) for x G S^. In case (i), 
2 — n < q < 2, we show that there is a p, —R < p < 0, such that p{x) > p^{x), for all x G S^. 
By a continuity argument, we show that the plane Xi = p can be moved continuously all 
the way to Xi = 0, and hence, p{x) increases as x approaches the origin from xi < 0. A 
similar argument can be made using planes xi = p > 0. Since the direction xi can be chosen 
arbitrarily we conclude that p is radially symmetric and decreasing about the origin. For 
case (ii), g > 2, a similar argument leads to p being radial and increasing about the origin. 
We now state and prove the result regarding the monotonicity of equilibria of ([I])-©- 

Theorem 3.2. Consider a bounded steady state p(x) of the aggregation model ([I])-® that 
is supported in a ball B{0,R) o/M". Then, p is radially symmetric and monotone about the 
origin. More specifically, we distinguish two cases: (i) 2 — n < q < 2, when p is decreasing 
about the origin, and (ii) q > 2, when p is increasing. 

Proof. We use the notations and symbols introduced in the preamble of the theorem. Cal- 
culate p(x) — p^(x), for X G S^. As p^{x) = outside B{0^, R), we have 

p(x) - p^(x) = p(x) > 0, for X G \ 5(0^ R). 

It remains to consider the case x G fl 5(0^, R). Denote by the complement of in 
5(0, i?), i.e., 

Sj = i?(0,i?)\S^, 

and calculate, using (!28|) . 

/^(^) = « / 1^ \\n-a Piy) dy + a -—1— -p(y) dy 

= / U_i|n-a ^(^) dy + a ^——p^{y) dy. 

In the above calculation we used \x^ — y^\ = \x — y\ to write the second integral in the 
right-hand-side as an integral with respect to y^ (subsequently relabelled y). 
Similarly, using ( 129|) . 



^/^(^) = " / Ir,^ -,An-a P(y) dy + a [ \^ _]An-a Pf^iy) 

Jt,^ f'^ y\ Js^nB{Qt^,R) F y\ 
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As is zero outside 5(0'^, i?), the second integrals in the right-hand-sides of the expressions 
for p and above can be extended to S^. Hence, we compute, for x G fl 5(0'^, i?), 

p(x) -p^{x) = aj^ (^—-L— - __L__^ (p(,^) _ p^(^)) rf^. (30) 

Case (i)2-n < g < 2 (or < a < n): For a; G S^ni?(0'', i?) and y G S^, < 
Hence, 

1 ^ - I ^ ^ In-. > 0' for all X G n S(0^ R) and y G S^. (31) 

Let us first assume that there exists po ^ 0) such that 

p{x) > PtMoi^), for all a; G S^o, (32) 

and prove that po can be extended all the way to the origin 0. Suppose by contradiction 
that po cannot be extended. Take p G {po, po + e), e > 0, and define 

K = e j:^ \ p{x) < pf,{x)}. 



For X G , using fl30|) and fl3T|) we find 



and hence, 

ilP/. - pIIl-(s-) < «IIPa< - pIIl-(s-) sup / ( , - , ^_\„_J (33) 

The function under the integral on the right- hand- side is integrable, as a > 0. Also, 
from fl30l) . we infer that p{x) > Pfj.{x) in the interior of S^j^, which implies that the clo- 
sure S~ of E~ has measure 0. As lim„ -.„„ S~ C S~ , we conclude that the measure of E~ 

Mo MO M^MO M MO ' M 

approaches as /i — )■ po- Therefore, we can choose e small enough such that 

a sup / -; ; -; ; dv < - . 

By (155]) we have Hp^t — p 11^00(2-) =0) which implies that S~ is empty, hence po is not maximal 
and we reached the desired contradiction. 

It remains to show that po G (— -R, 0) with the property fl32|) exists indeed. Note that fl32l) 
holds for po = — -R (see fl28l) and fl29|) ). An argument entirely similar to the one made in the 
previous step proves that the plane po = —R can be moved further to the right, while fl52|) 
still holds. This concludes the proof for case (i). 

Case (a) 2 < q (or a > n) follows similarly, the single difference being the sign of the 
integrand in f l5U]) . More precisely, instead of f l5T]) . the reversed inequality holds: 

^ ^ < 0, for all X G n 5(0'^, R) and y G S^. 



F — ?/|"~" |x^ — ?/| 
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A similar argument as that made for case (i) shows that p decreases as x approaches the 
origin from xi < 0. Hence conclude that in this case p is radially symmetric and increasing 
about the origin. □ 



Remarks 

a. Case q = 2 corresponds to constant solutions in a ball. These solutions were studied 
in detail in [2U] . 

b. Theoretical findings are in perfect agreement with numerical results — see Section |H 

c. The convexity of the steady state can be inferred easily from when g > 3 (the 
function is convex for a > 1). 

4 Numerical and asymptotic studies of radial equilib- 
ria 

We showed that a steady state p supported on a ball is necessarily radially symmetric 
and monotone. In this section we investigate further their properties using numerical and 
asymptotic methods. 

4.1 Numerical methods for the dynamic evolution and steady 
states 

In solving numerically the steady states ( 12^ or the dynamic evolution of solutions to ([I])- ([2]), 
the computational bottleneck is the evaluation of the integral operators. The methods we 
use are similar to those in [261 120] and are reviewed and extended below. 

Steady states. The equilibria are computed from ( 124]) by using the power method [T7] . 
whose convergence is guaranteed by Theorem 13.11 First, write the operator Tr given by (125|) 
in radial coordinates: 

TR-Pir) = ""r^. 1:^27 l\r'r'p{r')I{r. r')dr\ (34) 
sm" ^edO Jq 

where 

I(r, r') = (r^ + {r'f - 2rr' cos Oy^^"^ sin"-^ 0de. (35) 
Jo 

The steady states p are eigenfunctions of Tr corresponding to eigenvalue 1. Here the eigen- 
value problem consists in determining the eigenfunction p and the radius R of the support. 
The actual steady density is a constant multiple of this eigenfunction, where the constant is 
determined from the initial mass. 

By the scaling argument used to prove Theorem 13.11 , it is enough to find the spectral 
radius A of Ti, and its corresponding eigenfunction pi; the steady state p can then be 
calculated from ( 1271) and ( l26l) . 

Given an initial positive density p^^^ on [0,1], consider the iterative scheme (power 
method): 

p{m+i) ^ Tip^^Vll^ip^""^!!, = iitip^^^ii/iip^'")!!, (36) 
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where || ■ || can be any norm for functions on the unit interval [0, 1]. The sequence p*^'") 
converges to pi [17], and the spectral radius of Ti is given by the limit of X^"^\ as m goes to 
infinity. 

The integral operator in ( IMll is discretized by the trapezoidal rule. The computational 
complexity can be reduced to the calculation of the integration in r' only, where /(r, r') can 
be interpolated using the fact that I{r,r') = max(r, r')'?~^/i(s), where 

/•TT 

Ii{s)= {1 + -2scosey/^~^sm''-^ed9, s = min(r, r')/ max(r, r') G [0, 1], (37) 
Jo 

and Ji is computed at sample points on [0, 1]. 

The calculations do not present significant challenges, except for q E (2 — n, 3 — n], when 
/i(l) is unbounded, making the error in the trapezoidal rule uncontrollable. We managed 
to calculate this more singular case only in dimensions one and three, where J(r, r') can be 
obtained explicitly. Due to these computational difficulties, the singular hmit g — > 2 — n 
studied below by asymptotic methods is valid in all dimensions, but its numerical verification 
is done only in one and three dimensions. 

In three dimensions for instance, the kernel / from ( 135|) can be calculated explicitly: when 
r, r' 7^ 0, 



J(r, r') 



and when r' = or r = 0, 

/(r, 0) = 2r'?-2, l{0,r') = 2r"'-\ 

Using the explicit form ( l38l) . when q is in the singular range (—1, 0], and r ^ r', the integrand 
r''^I{r,r')p{r') in the integral operator is weakly singular. However, the part associated 
with this weak singularity can be approximated as the product of the weakly singular function 
|r — r'\~'^ and a smooth function. When the latter is approximated by linear or higher 
order interpolation, the whole integral can be calculated explicitly |18]. giving an accurate 
approximation of 

The steady states have been displayed in Figure [1] in dimension three, for a wide range of 
values of q. The results are perfectly consistent with the monotonicity properties stated and 
proved in Theorem 13.21 As noted from Figure [H equilibria display an interesting asymptotic 
behaviour as q ^ oo and g\2 — n. Asg— j-oo, the radii of the equilibria approach a 
constant, and mass aggregates toward the edge of the swarm. As g \ 2—n mass concentrates 
at the origin, as attraction becomes as strong as the Newtonian repulsion. We perform careful 
asymptotic studies in Section 14.21 and confirm and detail these observations. 

Dynamic evolution to equilibria. All numerical simulations we performed suggest that 
equilibria that solve f l24|) are global attractors for the aggregation model (II])- ([2]). To evolve 
dynamically the solutions to (II])-(l2]), we consider the model in characteristic form (see f|T6]l . 



{ {r+r')i-\r- 
qrr' 
Alnfi^ 
rr' \r—r' 
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(fT7|) . (ITT]) ). In radial coordinates the characteristic equations read 







/•oo /"TT 

/ (r)""V(^') / (r -r' cos 9){r'^ + r''^ -2rr' cos ey/^-^sin^'-^ededr' (39£ 

/•oo /"TT 

p-{n + q-2)ujn-i I (r')""V(0 / (r^ + r'^ - 2rr' cos^)''/^-^ sin"-^ ^d^dr' 



? = ^ f\rT~'p{r')dr' 

dp 

-17 = -P 

L -'o -'o 

(39b) 

The term associated with the singular repulsion (the Newtonian potential 0) in the right- 
hand-side of (I39ap is calculated by taking advantage of the fact that the corresponding kernel 
is the fundamental solution of the Laplace equation. 

Similar to the discretization of the integral from (IMll . the computational complexity in 
calculating the integrals from f l39ap and (I39bl) can be reduced by introducing the following 
auxiliary functions: Ji (defined by fl37j) ). 

/»7r 

l2{s)= / {l-scose){l + s'^-2scos9y/^-^sm''-^ed9, 
Jo 

/•TT 

Is{s)= / {s-cos9){l + s'^ -2scosey^'^-^sm''-'^ede. 

Jo 

Thus, the angular integrals in 6 in fl39ap and f l39bp become products of powers of r, r', and 
these auxiliary functions, with s = min(r, r')/ max(r, r'). Hence the double integrals in f l39ap 
and (139bp become single integrals in r' and are evaluated by trapezoidal rule. Due to the 
extra factor 1 — cos^ in the integrand, hi^) and 13(1) are bounded for any q > 2 — n, and 
can always be used in the trapezoidal rule. This observation reduces the total complexity in 
the computations to 0{N'^) per time step, where N is the number of spatial gridpoints in 
the radial coordinate r. Once the characteristic speeds in fl39l) are found, the equations are 
evolved in time by the classical fourth order Runge-Kutta method. 

Figures |2|^a) and Wljo) show simulation results in three dimensions, corresponding to 
g = 1.5 and q = 20, respectively. We plot the solution against the radial coordinate r. The 
initial data used in Figure |2] is 

p{x, 0) = (0.2 - 20|a;p + 1000|x|^) exp(-40|xn/c, (41) 

where c is a constant chosen to normalize the mass to one. The solutions approach as t — )■ 00 
the steady states studied in Section [3] and shown in Figure [H We note that for large q, the 
convergence near the origin tends to be slow. Based on numerical observations we conjecture 
that these equilibria are global attractors for solutions to ([I])- ([2]). In future work we plan to 
validate rigorously these observations regarding the global stability of the steady states. 

4.2 Asymptotic behaviour of equilibria as g — > 00 and q\2 — n 

In general, there is no explicit formula for the radius of the support R and the corresponding 
steady states p on B{0,R), governed by (^^. However, when q is large or close to the 
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(a) (b) 

Figure 2: Time evolution of a radially symmetric solution to the aggregation model ([I])-® 
with (a) q = 1.5 and (b) q = 20, in three dimensions, starting from initial data ( 14T]) . The 
solutions approach asymptotically the steady states studied in Section [3] (solutions to (l2ll) ) 
and shown in Figure [H Numerics with a variety of other initial conditions suggests that 
these equilibria are global attractors for the dynamics of ([I])-(l2]). 

singular limit 2 — n, the asymptotic behaviours can be obtained by perturbation expansions. 
To facilitate the exposition, we only consider the case R = 1 and solve the eigenvalue problem 

Tipi{x) = Api(x), for X e 5(0, 1), (42) 

where 

Tipi(x) = {n + q-2) [ \x - y\'^~^piiy)dy, (43) 

Jb(o,i) 

and pi vanishes outside 5(0,1). Provided we know the spectral radius A of Ti and its 
corresponding eigenfunction pi, the actual steady state p in can be recovered as p(r) = 
cpi{r/R), where R = A-V("+9-2) ^^^^ 

c is a constant related to the total mass (see the proof 
of Theorem 13. H in particular equations ( 126|) and (127|) ). 

In the rest of this section, the total mass for pi is assumed to be one. All asymptotic 
behaviours are investigated in one dimension first, to illustrate the essential techniques and 
features in a relatively simple setting, and then in higher dimensions. 

Asymptotic limit when g — j- oo. In this limit, the numerically computed eigenfunc- 
tions are observed to be concentrated near the boundary and the dominant contribution 
of the integral operator (H3|) comes from \y\ ^ 1. This motivates the following asymptotic 
construction. 

Start with the eigenvalue problem fl42p in one dimension (n = 1): 

Api(x) = (g - 1) y Jx - y\'^~'p,{y)dy. (44) 
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When X > 0, the above integral on the right hand side is dominated at ?/ ~ — 1, that is 

(g-l)y" \x-y\'^~'p,{y)dy^{q-l)p,{-l) j \x - yl^^^dy 

= pi(-l)((l+x)'?-i + (l-x)'^"i). (45) 

For X < we find the same expression, considering the symmetry pi(l) = p(— 1). Eval- 
uate (jS]) at a; = — 1 and use ( H5l) to approximate its right-hand-side. We find A ^ 2'?"^ 
and 

+ + (46) 

where the coefficient is chosen such that the approximated pi has total mass one. 

Approximations, as g — ?■ oo, to the actual steady states (with unit mass) and their 
support can be found from (l26l) and (1271) . Hence, R = A'^^'-'^^^^ ~ 0.5. However, this is not 
an accurate approximation, as it can be observed from Figure |3t^a), where we plot the radius 
R against q, as obtained from numerics (dots) and the above approximation (dashed line), 
the latter being referred to as the coarse approximation. To obtain the numerical results we 
used the methods described in Section 14.11 

The approximation of the eigenvalue A can be improved considerably if in the right- 
hand-side of ( l44l) one uses the expression of pi given by ( 146|) . Near y ^ —1, use pi{y) ~ 
q{l — yY~^/2'^'^^ to approximate the right-hand-side in (144|) . Then, evaluate at x = 1 to find 

The integral in the right-hand-side can be computed exactly and also, from fH^ . pi(l) ~ q/'i. 
We derive the refined approximation A ~ 2'^"^, with the corresponding radius, 

i?^2-(«-2)/(,-i)^ (47) 

The refined approximation, displayed as solid line in Figure |3l^a) shows an excellent agreement 
with the numerical results (dots). Note that R — )■ 0.5 (the coarse approximation), as g — )■ cxo, 
but the convergence is slow. 

Formally, the eigenfunction (H^ can be regarded as p^^^ obtained from the power method 
iteration (!36|) . by starting with a constant initial guess p^^^ = 1. The coarse approximation 
2'^"^ of the eigenvalue is exactly A*-^-*, while the refined approximation 2*^"^ is A*^^-*. Here the 
norm || ■ || is the function evaluation at x = 1, i.e., ||p|| = p(l). The same idea is applied 
below in higher dimensions, even though the expressions are much more complicted. This 
fact also illustrates the fast convergence of the iterative scheme fl36p for large q. 

We also compare the steady states, as obtained by numerics (see methods in Section 14. ip 
and asymptotics (expression fH6l) ). Figure |3t^b) shows an excellent agreement between the 
two solutions for g = 10 and g = 20. 

In higher dimensions, the eigenvalue problem can be approximated in a similar manner. 
Based on numerical observations, we assume that most contribution in the integral from fH5]) 
comes from the boundary. Hence, the eigenvalue problem for Ti is approximated by 

Xpi{x)^{n + q-2)pi{l) I \x -y\''-'^dy. (48) 
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(a) (b) 

Figure 3: Numerical and asymptotic solutions in dimension one, as g — ?■ oo. (a) Radius 
R of the steady states ( l24l) . The refined asymptotic approximation ( l47l) (solid line) agrees 
extremely well with the numerical solution (dots). The coarse approximation R = 0.5 
(dashed line) captures the (slow) asymptotic limit, as g — )■ oo. (b) Eigenfunctions pi of 
Ti computed by numerical (Section 14. ip and asymptotic (equation (146|1 ) methods. The 
agreement between the two sets of results is excellent. 



In general, the integral on the right-hand-side can not be integrated explicitly in the natural 
radial coordinates. However, when the origin is shifted to x, we have 



/ \x-yr'dy= [ Ivr'dy 

Jb[0,1) Jb{x,1) 



sin"-2 ede 

nu. 



sin'^-2^ / {r'Y+'i-^dr'de 

Jr2+r'2-2rr' cose<l 



+ g - 2) /g" sin"-2 OdO 



[n 



/ sin"-2 e{Vl - r2 sin^ ^ - r cos ey+'^-'^de. 
Jo 

(49) 



The approximation to the eigenvalue A is obtained by evaluating expression ( 14811 at r = 1, 
which gives 

= (n - l)c.„_i 2"+''-=^Beta ^^±|^ 

= (n - i)c.„_i2«+-^r(^)r(^i±pl)/r(n - 1 + 1), (50) 

where Beta and F are the beta and Gamma functions, respectively. 
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The radius of the support R = \ 2) ^^^^^ then be approximated as g — !■ 00 using 

Stirhng's approximation, 

T{z + 1) ~ V2^z'+^/^e-\ 

We refer to the outcome of this procedure as the coarse approximation and we plot the result 
in Figure Hl^a) (dashed line). Note the low accuracy of this approximation, when compared to 
the numerical calculation (dots). The coarse approximation approaches however, as g — )■ 00, 
the correct value 0.5, but the convergence is extremely slow. 
The eigenfunction pi is approximated from ( HHl) : 

Pi{r)^C sin"-2 e{Vl - r2 sin^ ^ - r cos ^)"+^-2(i^, (51) 
Jo 

where C is a normalization constant determined by the total mass. 

Similar to the ID case, using the explicit approximation f lHT]) of the eigenfunction, one 
can improve considerably the approximation of the corresponding eigenvalue A (and hence 
of R), by substituting (15T!) into (148|1 and evaluating at r = 1. The outcome of this improved 
procedure is referred to as the refined approximation and is plotted (solid line) in Figure |l](a). 
The agreement of the refined approximation with the numerical results is now excellent. 
There is also a very good agreement between the eigenfunctions computed numerically (see 
Section I^?T1) and their asymptotic approximation provided by (ISTl) — see Figure Ht^b). 

Remark. Since the integral expression fl5T]) for the approximation of pi is concentrated 
at 6* ^ TT, we can use Laplace's method to find the leading order when r is away from zero. 
Since pi(l) = C2"+«-3Beta (^, ^^±f^), and for 9 ^ n, 

- r cos ^ = (1 + r) [1 - r(7r - 9f/2 + ■■■], 
the eigenfunction pi can be further approximated away from the origin as 

pi(r) ^ C{1 + r)"+''-2 [ sin"-2 9(1 - r(7r - 9)^)"^''-^d9 

Jo 

POO 

= C(l + r)"+''"2 / Qn-2^-{n+g-2)re^^Q 

Jo 

= Cr~''"~"'"^/^(1 + r)^~^'^~'^. 

However, this expression breaks down for r near the origin and approximation (!5T]) is used 
instead in the plot from Figure Hl^b). 

To conclude, the asymptotic study as g — )■ 00 shows that the radii of support of the 
equilibria fl24p converge (slowly) to a fixed value R = 0.5, while the density concentrates 
to a (5-sphere of radius 0.5. We also point out that this asymptotic behaviour applies to all 
dimensions n. 

Asymptotic limit when g — )■ 2 — n. We write 

g = 2 - n + e. 
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Figure 4: Numerical and asymptotic solutions in dimension three, as q ^ oo. (a) Radius 
R of the steady states (!24|) . The refined asymptotic approximation (solid line) and the 
numerical solution (dots) agree extremely well. The dashed line represent the coarse approx- 
imation obtained from f l5Up and Stirling's formula. The radii R approach (very slowly) the 
constant 0.5, as g — )■ oo. (b) Eigenfunctions pi of Ti computed by numerical and asymptotic 
(approximation f l5T|) ) methods. There is very good agreement between the two solutions. 



and perform an asymptotic study in the small e regime. We start with solutions in one 
dimension first, where the eigenvalue problem reads 



XePlix) 



1 

el \x- y\'~^pl{y)dy, 



q-1. 



(52) 



Here the subscript and superscript e is used to emphasize the dependence of the eigenvalue 
A and the corresponding eigenfunction pi on e. The asymptotic expansion suggested from 
numerical simulation is 

A, = Ao + Aie + Aae^ + ■ • • , (53a) 

pl (x) = p(°) (x) + ep(^) (x) + e2p(2) (x) + ■ ■ • . (53b) 

The kernel \x — is not integrable in the limit when e — 0, and thus we can not 
substitute the formal expansions into f l52|) directly. Hence, it is not possible to carry out a 
straightforward expansion. However, using the fact that 



\x — '^dy = (1 + xY + (1 — x) 



the governing equation ( 152|) can be written as 
[A,-(l-x)^-(l + x)^]p^i(x) = 



el \x — y 



[Pl{y)-Pl{x)]dy. 



(54) 
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Now using the asymptotic expansions fl53ap and (]53bp . we obtain 



0(1) : Ao - 2 = (55a) 



0(e): j \y-x\~\p^'>\y)-p^'>\x))dy-{\i-\n{l-x^))p^''\x)) = (55b) 
0{e') : \y - x^^p'^'^y)) - p^'\x)))dy - (Ai - ln(l - x'))p^'\x)) 



The first equation fl55ap yields Ao = 2 and the second equation (]55bp is an eigenvalue problem 
for the limiting profile p^^\ Integrating the equation ( I55bl) with respect to x, we can get an 
alternative expression for the eigenvalue Ai as 

= J-^i" ^ ^ ^ (56) 



Even though we do not expect any explicit solutions to f l55bp . we can discretize and solve 



it using inverse iteration [22]. The condition p*^'^^(±l) = 0, motivated from the numerical 
calculation of the steady states, is used to get rid of the singularity on the boundary. The 
initial data for the inverse iteration can be taken as the steady state p\ calculated numerically 
from fl52|) with e close to zero (and initial guess of the eigenvalue from f l56|) with p(°^ replaced 
by p\). This inverse power iteration normally converges in just a few steps. 

The solution p^^^ of fl55cp is more challenging. Since this first order correction p^^) does 
not provide much insight into the problem, we focus instead on the second order correction 
A2 of the eigenvalue A. By the solvability condition for f l55cp . the right hand side of f l55cp is 
orthogonal to p'-^^ giving 



In 


y-x\ 




y-x 





X2 I [p^''\x)] dx= I _ Si (P^^Hy) - P^^Hx))p^^\x)dydx 

[p(°)(x)]'rfa;. (57) 



ln^(l -x) +ln^(l + x) r_tm, .i2 



Figure [5]^a) shows the normalized steady states pi for e = 1,0.5,0.2 (corresponding to 
q = 2, 1.5, 1.2), as computed numerically from the eigenvalue problem ( l52l) using the methods 
described in Section 14.11 In the same figure we also plot (plain solid line) the leading order 
term p*-°^ of the expansion (I53bp . obtained by solving ( I55bp with the inverse iteration method. 
The plot confirms that the equilibria pi approach the limiting profile p*-*^^ as e — )■ 0. 

In Figure [5](b) we plot (dots) the eigenvalues A^ for various values of e G (0, 1] (corre- 
sponding to g G (1,2]), computed directly from ( l52l) using the power method (Section l4.ip . 
The dashed line represents the asymptotic approximation of A^ from f l53al) that includes sec- 
ond order corrections (Ai is computed from ( l56l) and A2 from ( 1571) . The agreement between 
the two sets of results for small e's is excellent. 

At a closer inspection, it becomes clear that the expansion (]53bp may be non-uniform 
near the boundary x = 1. In Figures [6](a)-(b) we plot pi evaluated at x = and x = 1, 
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(a) (b) 

Figure 5: (a) Normalized steady states p\ for e = 1, 0.5, 0.2 (g = 2, 1.5, 1.2) computed directly 
from (152!) using the numerical methods from Section l4m The plain solid line represents the 
leading order term p^"-* of the expansion (]53bl) . Equilibria pi approach the limiting profile 
p^^^ as e — ?■ 0. (b) Eigenvalues for various values of e G (0, 1] (g G (1, 2]), computed directly 
from f l52|) (dots). The dashed line is the asymptotic approximation of from fl53ap valid at 
order O(e^). Note the excellent agreement between numerics and asymptotics for small e [q 
close to 1). 



respectively, for different values of e. The eigenf unctions pi are computed directly from 
— see also Figure [5l^a). At the origin (Figure |6]|^a)) we find a linear dependence on e, as 
the approximation pl{0) = p^^\0) + ep^^\0) + ■ ■ ■ is uniform. However, at x = 1, we find 
pl{l) ~ V0.253e + O(e2), hence the expansion fl53bl) is non-uniform near the boundary 
(Figure Et^b)). To find a valid asymptotic expansion near a; = 1 one has to introduce a 
boundary layer and perhaps use the method of matched asymptotics to relate the inner and 
outer expansions. We do not pursue this direction here. 

The situation in higher dimensions is similar. We write the eigenvalue problem (142!) as 

(A, - h{r))pl{r) = e I \x - yr^'iPliy) - pl{x))dy, r = (58) 

Jb{0,1) 

where we used again subscripts and superscripts to emphasize the dependence on e = q+n—2, 
and the auxiliary function is defined by 

k,{\x\) = e / Ix-yl'-'^dy. 
Jb{o,i) 

Equation (158|1 is the analogue of (154|) from ID. 

The auxiliary function k^{r) can be simplified by shifting the origin (see the calculation 
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Figure 6: (a) Validation of the asymptotic expansion (]53bp away from the boundary. At 
X = 0, p\{0) has a hnear dependence on e of the form 0.5944 — 0.0936e, as the expansion (]53bp 
is uniform near the origin, (b) At x = 1, has a square-root dependence on e of the 

form V0.2530e. There is a boundary layer near x = 1 where the expansion (I53bl) is no longer 
vahd. 



leading to (l49l)). with the result: 



sin'^-2 Ode 



nUr. 



sin"-2 Ode 



sin'^-^ eiVl-r^ sin' 9 - r cos eYdO 



sin'^-^ e(l + e ln( Vl - sin^ ^ - r cos 6) 



ln2( Vl - r2 sin^ ^ - r cos ^) + 0(e^)) ci^, 



where we also used a Taylor expansion in e to expand (y 1 — sin^ ^ — r cos^)^. The 0(e) 
term can be simplified using the the following calculation: 

sin"-2 e ln(\/l - r2sin2^ - r cos 



7r/2 



sin"-2 ^ 



ln( V 1 — ''"^ sin^ 6 — r cos ^) + ln( y 1 — sin^ 6* + r cos 



de 



ln(l-r2) / sin"-^^^^. 



Hence, 



K{r) =nun + '^ ln(l - r')e + k^'\r)e^ + 0(e=^), 



(59) 
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where 

A;(2)(r) = ^ r sin"-2 ^ In^r - r2 sin^ ^ - r cos e)de. 

2 sin'' ^rf^ Jo 

Since the integrand in (l58l) is integrable provided is Holder continuous, by introducing 
the formal asymptotic expansion (|53ll in (|58l) and using (|59l) . one finds 

0(1) : Ao — nUn = 

0(6): / |^_^|-(p(0)(3;)_p(0)(^))rf^ 

JB(0,1) 

0(e2): / |^_a;|-"(pW(t/)-p«(x))dy 

JB{0,1) 

(A2-A;(2)(|a:|))p(o)(x)- / ^ 

Jb{0,1) 

Therefore, the leading order of the eigenvalue is Aq = ncUn and the limiting steady state 
p^'^^ and the first order correction Ai can be solved by inverse iteration. The second order 
correction A2 can be obtained from the solvability condition: 

A2 /" [p^''\x)\^ dx - f k^''\\x\)[p^^\x)f dx 
= [ [ ly-xl-^'lnly-xllp^^^y)- p^^\x))]p^°\x)dydx (61) 

In Figure [7] we test the asymptotic results in three dimensions. The qualitative behaviour 
of the steady states p\ and their corresponding eigenvalues A^ is similar to what has been 
observed in one dimension (Figure E]). 

Figure [Tl^a) shows the normalized steady states p\ for e = 2, 1,0.5,0.2 (corresponding to 
q = 1,0,-0.5,-0.8) obtained by solving numerically (142|) . along with the limiting profile 
p*^*^-* (plain solid line) found by inverse iteration from (160bp . The asymptotic results are 
confirmed, as equilibria pi approach the limiting profile p'-^^ for e — >■ 0. In Figure W(Jo) 
we test the asymptotic expansion f l53al) at order O(e^) (dashed line) against the numerical 
solution of (H2I) (dots). The agreement is excellent for small e or equivalently, for q close to 
the critical value —1. 

Remark. An implication of the above study is that, in all dimensions, the radius R = 
^-i/{n+q~2) _ yauishes exponentially fast as g \ 2 — n, and as a result, the true steady 
state p converges to a Dirac delta function. 

5 Discussion 

We have studied the aggregation model ([1]) with potentials K that contain short-range 
Newtonian repulsion and long-range power-law attraction. The main merit of the family 
of potentials considered here is that it leads to solutions which have biologically relevant 



(60a) 

-(Ai-^na;„ln(l-|x|2))p(°)(x) = 

(60b) 

- (^Ai - ^ncOnln{l - \x\^)^ p^'\x) = 

- In \y - x|(p(°ni/) - P^'^\x))dy. 

(60c) 
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Figure 7: (a) Equlibria p\ in three dimensions for e = 2,1,0.5,0.2 (g = 1,0,-0.5,-0.8) 
computed numerically from P2l) (see Section HTTl) . The plain solid line is the limiting profile 
p^^^ found by inverse iteration from fl60bl) . As e — j- 0, p\ approaches the limiting profile, 
confirming the asymptotic expansion, (b) Comparison between numerics and asymptotics in 
three dimensions: eigenvalues obtained numerically from ( 142|) (dots) and from the asymptotic 
expansion f l53ap valid at order O(e^). The two sets of results agree very well for small e (g 
close to —1). 



features, such as finite densities, sharp boundaries and long lifetimes [22] • Finding such 
solutions to model ([T]) has been indicated as a "challenge" in previous works [531 121] and the 
literature addressing this issue has been very scarce. 

Well-posedness of solutions to ([1]) was studied by analogy with incompressible fluid equa- 
tions [30], using the Lagrangian (particle) formulation of the model. The main object of 
the present work, i.e., equilibria supported on a ball, was investigated through a variety of 
analytical, numerical and asymptotic methods applied to the integral equation 0241) . We 
derived existence and uniqueness of such equilibria using the Krein-Rutman theorem and 
established qualitative properties such as monotonicity and radial symmetry, by the method 
of moving planes. 

The numerical results confirm the analytical findings and also suggest that the equilibria 
studied in this work are global attractors for the dynamics of ([T]). We formulate this obser- 
vation as a conjecture and we plan to address it in future work. A possible approach is to 
use the fact that the aggregation equation ([1]) represents a gradient flow with respect to the 
energy 

E[p] = 7; / -y)p{x)p{y)dydx. 

However the energy is not convex, therefore its global minimizers cannot be characterized 
easily. Some recent progress in this direction was done in [H] in the context of aggregation 
models with long-range attraction and quadratic diffusion. 

The asymptotic results revealed some very interesting features of solutions to In 
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particular, as the exponent q of the power-law attraction approaches oo, the radii R of the 
support approach a constant value 0.5, and the density concentrates on a (5-sphere. Distri- 
butions on spheres (uniform, as well as surprisingly complex patterns) have been studied 
recently [371 [1] using potentials with power-law repulsion and attraction. Concentrations on 
5-spheres typically represent equilibrium solutions of the aggregation model. It is not sur- 
prising in fact that the non-convex energy E has multiple stationary points, adding to the 
difficulties in studying its equilibria, as indicated before. Aggregations on spherical shells 
could be stable or unstable, depending on the exponents of the repulsive and attractive 
power-laws [271 [1]. Choosing the repulsion component in Newtonian form, as in this pa- 
per, seems to rule out concentrations on spheres from the possible asymptotic behaviours of 
solutions to ([1]), except in the limit g — oo. 

The limit g — )■ 2 — n of solutions to ( 124]) is interesting for its own sake. Weakly singu- 
lar integral operators are subjects of many articles and textbooks (see [SB] and references 
therein), but a careful asymptotic study of the eigenvalue problem ( 121]) . as the singularity 
approaches the critical value 2 — n, is missing from the literature. We studied the scaled 
problem f H2|) for R = 1 and showed that eigenvalues A approach a constant, while the cor- 
responding eigenfunctions approach a limiting profile p*-"-*, as g — 2 — (see Figures [5] and 
[7]). Consequently, from the scaling ( l26l) -(l271). solutions to ( l24|) approach a Dirac 5 in the 
limit. The findings are consistent with works that consider blow-up in aggregation models 
with purely attractive potentials [T9l H] , in particular recent works that consider Newtonian 
potentials [7]. 

Finally, we want to comment briefiy on the biologically unrealistic feature of the potential 
([2]), that is, the growth of attraction with distance, when g > 0. As discussed in more detail 
in [2U], the dynamics of ([T]) remains unchanged if the potential is modified in an arbitrary 
way outside a ball of M", with a sufficiently large radius that depends on the initial conditions 
only. This can be inferred from the property of the density to have uniform (in time) support; 
this property was shown to hold for g > 1 and it is believed to hold for all g > 2 — in fact. 
Provided the radius of the support of the density p{x, t) is bounded by -Rmax, where i?max 
depends only on the initial conditions, but not on time, the potential K{r) can be taken to 
be zero (or exponentially decaying) for r > 2i?jnax; without changing the dynamics. We refer 
the reader to [20] for a numerical illustration of this issue. 

6 Appendix 

Proof of Theorem 2.1 We show that the operator J-" : Ol — ?■ i3 is bounded, i.e., 
||J-'(X)||i ,y < oo, for all X G Ol- We outline the main steps and refer to Chapter 4 [3U] for 
details and various technical calculus inequalities. 
Write F{X) as 

F{X) = voX. 
Using the calculus inequality (Lemma 4.1 [30] ) 

— II^IIl^^I^It ~i~ I^l7 II ^ II ; (62) 

we estimate: 

\\FiX)\\^^^ < Mlo. + (WVvWloo + \\/v\^) ||X||i,^. 
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As < L for X G Oli it remains to bound ||t>||ioo, ||Vt'||Loo, and iVt"]^. 

We first inspect the repulsion component of v (see ( l6b1) . ( !6cl) and ( fT9l) ) and estimate 
II A; * pIIoo, ||-P[p]||l°° and |P[p]|^, where P = VA; and P[p\ is the principal-value SIO: 

P[p](x)=PV y P(a: - 

The first of these terms can be bounded as follows: 

I [ 1 I [ 1 

nUn J\x-y\<l F - y\ J \x-y\>l F " 1/1 

„ „ M 

< \\p\\l°° H ■ (63) 

ncun 

To bound P|p] we follow the proof of Lemma 4.6 [30] • We split the integral: 

|P[p](a:)| =PV [ P{x-y)p{y)dy+ f P{x - y)p{y)dy . 

J \x~y\<e J\x—y\>e 



h{x) hix) 

The kernel P has mean- value zero on the unit sphere, 

Pds = 0, 

x\=l 

which enables us to rewrite /i as 

Ji(x) = PV / P{y){p{x -y)- p{x))dy. 
■J\y\<^ 

Hence 



|/i(a;)| < / P{y) r\- wydy 



\y\<e \y 



7 



< c|p|, / lyr-'-'dy, 

and 

\h{x)\ < c|p|^e^. 

Here and below, c denotes a generic constant. For the second integral, we estimate 
|/2(a;)| < / \P{x - y)\p{y)dy + / \P{x - y)\p{y)dy 

•J e<\x—y\<l J\x—y\>l 

< c||p||i-log 0^ +M. 

Conclude 

||P[p]||l=o <c(|p|,e^ + ||p|Uoolog(l/e)) + M. (64) 
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Finally, using the same argument as in Lemma 4.6 [30], one can show 

\P[p\\,<c\p\,. (65) 

Estimate the seminorm |p|^ as in the proof of Proposition 4.2 [30], by using ( ITOl) and calculus 
inequality f l62|) : 

|pl7 < IIpo o^^li-ldet V^X^^I^ + Ipo oX~^|^||det V^.X"^||ioo. (66) 
Using the calculus inequalities 

|/oF|,< |/|,||VF||Ioo, (67) 
\\Y~%,,<c\\Y\\l^:-\ (68) 



from Lemmas 4.3 and 4.2 [30], respectively, estimate 

< IP0I7 ||V, 

< c|po|7l|-^lll7 



P0O^"'l7<IP0UI|V.X-l"^ 



7 _ IFOI7 II * X^*- II 

7(2n-l) 



and 

|det V.X^^I^ < c||X||f7\ 
As X G Cl, ||det Vx-X-1||loo < L and < L. Return to ([66]) to find: 

\pIi < C{L) (IIpoIU- + IP0I7) • 

Using the above bound on |p|^ and the uniform bound Oil W p\ \ r.oc derived in Section 12. 1^ 
we infer from flB5]) - fl5S]) that the repulsion component of J-" yields a bounded operator. The 
singularity of the attraction component is milder than that of the repulsion and does not 
break the existing estimates for the singular repulsion kernel k. We conclude that J-" is 
bounded. 

To prove that J-" is Lipschitz continuous, we show that J^'{X) is bounded as a linear 
operator from Ol to S, i.e., ||J-''(X)|| < 00, for all X G 0^. Calculate J^'{X) using (fT7|) and 
3: 

:F\X)Y = ^F{X + eY) 

Vf{X{a) - X(/3))(F(«) - Y{f3))po{f3)df3, (69) 

where we suppressed the time dependence for convenience. 

To estimate the component of J^'{X)Y due to repulsion one could follow the proof of 
Lemma 4.10 in [SU]. The attraction component is milder and does not break the estimates. 
It can be shown that 

\\r{X)Y\U^, < C{L) (IIpoIIl- + IP0I7) nii,7> 
which proves the boundedness in the operator norm of J^'{X). 
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A key observation is that the term Y{a) — Y{I3) in (169!) compensates for the singularity 
of V/. We present here the estimate of || ||ioo, which requires in fact some changes 
to the proof of Lemma 4.10 [30]. We do not present the estimates of || Vq,J-''(X)F||loo and 
^, we refer instead to the calculations in |30j. 

Use a change of variable, x = X{a), y = X{P), and split the repulsion component of 
riX)Y from §^ 

Vk{x - y){Y{X-\x)) - Y{X-\y)))p{y)dy = 



x~y\<l J \x—y\>l 
Jl J2 

The procedure is similar to what we did to derive (p5]) . From mean- value theorem, we have 

\Y{X~\x)) - Y{X-\y))\ < \\{VaY o X-')V,X-'\\lo.\x - y\. 
As \Wk{x)\ < clxl-"", 

\Ji\<c\\{VaYox~')v.x~'\\L^\\p\\L^ [ - — K^dy 

J\x-y\<i F y\ 

<C(L)||V,F|U=o, 

where we also used (168|) . X & Ol, and the uniform bound on ||p||loo in the second inequality. 
The outer integral satisfies 

\J2\ < c\\{V^Y o X-')V,X-'\\l^ [ p{y)dy 

J\x-y\>l 

<C(L)M||V„r||Loo, 

where we used conservation of mass. 

Combine the two estimates for Ji and J2 and argue that the attraction component would 
not break these estimates, to find 



\riX)Y\\L^<CiL)\\Y\\,, 
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